/* 
 * File:   PreProcess.cpp
 * Author: song
 * 
 * Created on 2010年10月28日, 下午12:26
 */

#include "PreProcess.h"

extern int hdims;
extern int minT;
extern int numFrame;

PreProcess::PreProcess(Camera *camera)
{
    this->alpha = 0.03;
    this->cam = camera;
}

PreProcess::PreProcess(const PreProcess& orig)
{
}

PreProcess::~PreProcess()
{
}

//对外接口，处理图像
void PreProcess::Process(IplImage *image, Camera *camera)
{
    this->cam = camera;

    //图像灰度化
    this->ToGray(image, cam->pGray);
    //暂时使用高斯滤波
    this->ToSmooth(cam->pGrayMat, cam->pGrayMat, 1);

    if(numFrame == 1){
        cvZero(cam->pBkMat);
        //cvCopy(pGrayMat, pBkMat);
        //cvCopy(pGrayMat, pFrMat);
    }

    //非处于背景初始阶段
    if(numFrame > 50) {
        //帧差 二值化得到前景图
        this->ToBinary(cam->pGrayMat, cam->pBkMat, cam->pFrMat);
        
    }
    
    this->UpdateBK();
    //第五十一帧的时候前景矩阵值不为0和255原因不明
    if(numFrame > 50){
        //cvShowImage("bk", pFrImg);
        
        this->ToOpen(cam->pFrImg, 3);
        //this->ToClose(pFrImg, 5);
        //cvShowImage("bk", pFrImg);
//        if(numFrame > 51){
//        CvMat *mat = cam->pFrMat;
//        for(int y = 0; y < mat->height; ++y)
//        {
//            for(int x = 0; x < mat->width; ++x)
//            {
//                cout<< *(mat->data.fl + y * mat->width + x) << "  ";
//            }
//        }
//        int i  = 1;
//        }
    }
}
/////////////////////////////////////////////////////////////////////////////

//灰度化图像至前景中，并生成对应的矩阵
void PreProcess::ToGray(IplImage* image, IplImage* des)
{
    cvCvtColor(image, des, CV_BGR2GRAY);
    cvConvert(des, cam->pGrayMat);
}

///////////////////////////////////////////////////////////////////
void GuassFilter(CvMat *src, CvMat *des)
{
    cvSmooth(src, des, CV_GAUSSIAN,3,0,0);
}

void MeanFilter(CvMat *src, CvMat *des)
{
    int s = 5, r = s >> 1;
    double sum;
    if(src->type != CV_32FC1) return;

    try
    {
        for(int i = r; i < src->height - r; ++i)
        {
            for(int j = r; j < src->width - r; ++j)
            {
                if(j == r)
                {
                    sum = 0;
                    for(int x = -r; x <= r; ++x)
                    {
                        for(int y = -r; y <= r; ++y)
                        {
                            sum += cvmGet(src, i + x, j + y);
                        }
                    }
                }
                else
                {
                    for(int y = -r; y <= r; ++y)
                    {
                        sum += cvmGet(src, i + r, j + y) - cvmGet(src, i - r - 1, j + y);
                    }
                }
                cvmSet(des, i, j, sum / s / s);
            }
        }
    }
    catch(...)
    {
        printf("均值滤波失败！");
    }
}

//type为1使用高斯滤波，type为2使用均值滤波,内存泄漏
void PreProcess::ToSmooth(CvMat *src, CvMat *des, int type)
{
    switch(type)
    {
        case 1:
            GuassFilter(src, des);
            break;

        case 2:
            CvMat *tmp = cvCreateMat(src->rows, src->cols, src->type);
            MeanFilter(src, tmp);
            cvCopy(des, tmp);
        break;
    }

}

////////////////////////////////////////////////////////////////////
//根据灰度直方图统计计算适应阈值,实验结果阈值计算偏大，原因不明，代码实现可能有问题
//int OSTU1(CvHistogram *hist, Camera *cam)
//{
//    //转换为hsv图片格式
//    cvCvtColor(cam->pFrame, cam->hsv, CV_BGR2HSV);
//
//    //提取亮度分量
//    cvSplit(cam->hsv,cam->hue, 0, 0, 0);
//    //计算灰度直方图
//    cvCalcHist(&cam->hue, cam->hist);
//
//    float histArr[256];
//
//    //hdims = 256, 统计每个灰度级别下的像素个数
//    for(int i = 0; i < hdims; i++){
//        histArr[i] = cvQueryHistValue_1D(hist, i);
//    }
//    histArr[0] = 1;
//
//    float u0, u1, w0, w1;
//    int count0;
//    //t为阈值 maxT为方差最大时的阈值
//    int t, maxT;
//    //方差和最大方差
//    float devi, maxDevi = 0;
//    int i, sum = 0;
//
//    for(i = 0; i < hdims; i++)
//        sum += histArr[i];
//    for(t = 0; t < 255; t++){
//        u0 = 0;
//        count0 = 0;
//        for(i = 0; i <= t; i++){
//            u0 += i*histArr[i];
//            count0 += histArr[i];
//        }
//        u0 /= count0;
//        w0 = (float)count0/sum;
//
//        u1 = 0;
//        for(i = t+1; i < hdims; i++)
//            u1 += i*histArr[i];
//
//        u1 /= sum - count0;
//        w1 = 1-w0;
//
//        devi = w0*w1*(u1-u0)*(u1-u0);
//        if(devi > maxDevi){
//            maxDevi = devi;
//            maxT = t;
//        }
//    }
//    if(maxT < minT) maxT = minT;
//    //delete(histArr);
//    return maxT;
//}

//初步实验效果还不错
int OSTU2(IplImage *image) {
    int w = image->width;
    int h = image->height;

    unsigned char *np; // 图像指针
    unsigned char pixel;
    int thresholdValue = 1; // 阈值
    int ihist[256]; // 图像直方图，256个点

    int i, j, k; // various counters
    int n, n1, n2, gmin, gmax;
    double m1, m2, sum, csum, fmax, sb;

    // 对直方图置零...
    memset(ihist, 0, sizeof (ihist));

    gmin = 255;
    gmax = 0;
    // 生成直方图
    for (i = 0; i < h; i++) {
        np = (unsigned char*) (image->imageData + image->widthStep * i);
        for (j = 0; j < w; j++) {
            pixel = np[j];
            ihist[ pixel]++;
            if (pixel > gmax) gmax = pixel;
            if (pixel < gmin) gmin = pixel;
        }
    }

    // set up everything
    sum = csum = 0.0;
    n = 0;


    for (k = 0; k <= 255; k++) {
        sum += k * ihist[k]; /* x*f(x) 质量矩*/
        n += ihist[k]; /* f(x) 质量 */
    }

    if (!n) {
        // if n has no value, there is problems...
        //fprintf (stderr, "NOT NORMAL thresholdValue = 160\n");
        cout<<"NOT NORMAL thresholdValue"<<endl;
        thresholdValue = 80;
    }

    // do the otsu global thresholding method
    fmax = -1.0;
    n1 = 0;
    for (k = 0; k < 255; k++) {
        n1 += ihist[k];
        if (!n1) {
            continue;
        }
        n2 = n - n1;
        if (n2 == 0) {
            break;
        }
        csum += k * ihist[k];
        m1 = csum / n1;
        m2 = (sum - csum) / n2;
        sb = n1 * n2 * (m1 - m2) * (m1 - m2);
        /* bbg: note: can be optimized. */
        if (sb > fmax) {
            fmax = sb;
            thresholdValue = k;
        }
    }

    if(thresholdValue < 40) thresholdValue = 40;
    return thresholdValue;
}
void PreProcess::ToBinary(CvMat *frame, CvMat *bk, CvMat *fr)
{
    cvAbsDiff(frame, bk, fr);

    //查看矩阵
//    float *a = fr->data.fl;
//    for(int i = 0; i < fr->rows; i++){
//        for(int j = 0; j < fr->cols; j++){
//            cout<<a[i*fr->cols+j]<<" ";
//        }
//        cout<<endl;
//    }
   
    //获取阈值
    int t = OSTU2(cam->pGray);
    t /= 2;
    if(t < 40) t = 40;
    cvThreshold(cam->pFrMat, cam->pFrImg, t, 255.0, CV_THRESH_BINARY);
    cvConvert(cam->pFrImg, cam->pFrMat);
//    if(numFrame > 50){
//        CvMat *mat = cam->pFrMat;
//        for(int y = 0; y < mat->height; ++y)
//        {
//            for(int x = 0; x < mat->width; ++x)
//            {
//                cout<< *(mat->data.fl + y * mat->width + x) << "  ";
//            }
//        }
//    }
//    cvConvert(cam->pFrMat, cam->pBkImg);
//    //自适应二值化
//    cvAdaptiveThreshold(cam->pBkImg, cam->pFrImg,255.0,CV_ADAPTIVE_THRESH_MEAN_C,
//                        CV_THRESH_BINARY);
    
//    CvScalar val0 = cvScalar(0), val1 = cvScalar(255);
//    for(int i = 0; i < cam->pFrImg->height; ++i)
//    {
//        for(int j = 0; j < cam->pFrImg->width; ++j)
//        {
//            if(cvGet2D(cam->pFrImg, i, j).val[0] != 0)
//            {
//                cvSet2D(cam->pFrImg, i, j, val0);
//            } else {
//                cvSet2D(cam->pFrImg, i, j, val1);
//            }
//        }
//    }
    //cvNot(cam->pFrImg, cam->pFrImg);

}

void PreProcess::ToOpen(IplImage* image, int num)
{
    //开操作
    cvErode(image, image, 0, num);
    cvDilate(image, image, 0, num);
}

void PreProcess::ToClose(IplImage* image, int num)
{
    //关操作
    cvDilate(image, image, 0, num);
    cvErode(image, image, 0, num);
}

///////////////////////////////////////////////////////////////////////
//获得去除当前目标物体像素后的tmpBk

//输入当前帧灰度图像和目标图片,效果不好放弃使用
//void GetTmpBk(IplImage *image, CvMat *tmpBk)
//{
//    //使用wrapper方便像素访问
//
//    BwImage myGray;
//    ///本机实验摄像头返回的图像depth为IPL_DEPTH_8U，这里默认是这样吧
//    try{
//    myGray = image;
//    }
//    catch(...){
//        cout<<"背景更新GetTmpBk函数中，图像格式出错"<<endl;
//        return ;
//    }
//
//    float *fg = pFrMat->data.fl;
//    float *bk = pBkMat->data.fl;
//    float *tmp = tmpBk->data.fl;
//    int n = pFrMat->cols;
//
//    for(int i = 0; i < pFrMat->rows; i++) {
//        for(int j = 0; j < pFrMat->cols; j++){
//            //cout<<fg[i*n+j]<<"  ";
//            //为目标区域时用之前的背景像素,255为白
//            if(fg[i*n+j] >= 200) {
//                tmp[i*n+j] = bk[i*n+j];
//            }
//            //为背景区域时用当前帧的像素
//            else{
//                tmp[i*n+j] = myGray[i][j];
//            }
//        }
//    }
//
//}
void PreProcess::UpdateBK()
{
    if(numFrame <= 50) {
        cvAddWeighted(cam->pBkMat, 1, cam->pGrayMat, 0.02, 0, cam->pBkMat);
    }
    else if(numFrame > 50){

        //获取基于Surendra的背景-------------效果不好，一旦运动物体被背景录入，而物体一旦离开原来的位置，则会被永远检测成前景
        //GetTmpBk(pGray, tmpBkMat);
        ////////调试代码////////////////
        //cvConvert(tmpBkMat, tmpBk);
        //cvShowImage("tmpBk", tmpBk);
        //////////////////////////////

        //背景更新
        //cvRunningAvg(tmpBk, pBkMat, 0.05, 0);
        //cvRunningAvg(cam->pGrayMat, cam->pBkMat, 0.03, 0);
        //if(numFrame % 2 == 0)
        //cvCopy(cam->pGrayMat,cam->pBkMat);
    }
    
    //帧计数
    numFrame++;
}